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ABSTRACT 

The  eigenvector  solution  of  the  time- invariant  matrix  Riccati  equation 
is  discussed.   The  coefficient  matrix  of  the  canonical  equation  is  allowed 
to  have  multiple  eigenvalues,  namely,  the  matrix  could  be  either  derog- 
atory or  defective.   The  solution  of  matrix  Riccati  equation  is  then  cal- 
culated from  a  part  of  similarity  transformation  which  should  reduce 
the  coefficient  matrix  to  the  Jordan  canonical  form. 


Ill 
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1 .   INTRODUCTION 

It  is  well  known  [1,2]  that  the  linear  regulator  problem  with  quadratic 
cost  functional  is  reduced  to  the  problem  of  the  matrix  Riccati  equation. 

Although  the  solution  of  the  matrix  Riccati  equation  exists  and  is 
unique  [2,3],  it  is  not  easy  to  obtain  a  numerical  solution  for  the 
high-order  system. 

One  of  the  numerical  methods  is  the  eigenvector  solution  which  has  been 
studied  by  several  authors  [U,5,6].   In  their  studies,  however,  the  coeffi- 
cient matrix  of  the  canonical  equation  is  assumed  to  have  distinct  eigen- 
values, which  is  a  restriction  to  the  method  of  eigenvector  solutions. 
Recently,  Martensson  [7]  discussed  the  case  in  which  the  coefficient  matrix 
is  nondiagonable.   Similarly,  in  this  paper,  the  assumption  of  distinct 
eigenvalues  will  also  be  removed  so  that  the  coefficient  matrix  is  allowed 
to  have  the  multiple  eigenvalues,  namely,  the  matrix  is  either  derogatory 
or  defective.   To  do  this,  we  explicitly  construct  the  similarity  transfor- 
mation which  should  reduce  the  coefficient  matrix  to  the  Jordan  canonical 
form,  either  diagonal  or  nondiagonal,  and  then  the  solution  of  the  matrix 
Riccati  equation  is  calculated  from  a  part  of  the  similarity  transformation 
matrix. 

Since  matrix  Riccati  equations  of  fairly  large  sizes  appear  in  many 
engineering  applications  and  since  the  solution  algorithm  can  be  easily 
constructed  such  that  it  becomes  suitable  for  a  parallel  computer  [8],  we 
provide  in  Appendix  II  the  basic  codes  written  in  ILLIAC  IV  language, 
GLYENIR. 


2.   FORMULATION  OF  THE  PROBLEM 

In  this  section,  the  source  of  the  matrix  Riccati  equation  and  its 
relation  to  optimal  control  problems  are  summarized  in  short  for  our 
further  discussions  [1,2,3] . 

Let  us  consider  the  linear  time- invariant  dynamical  system 

x(t)  =  Ax(t)  +  Bu(t) 

y(t)  =  Cx(t)  (2.1) 

x(tQ)=  xQ 
where  A,  B,  and  C  are  n  x  n,  n  x  r,  and  m  x  n  constant  matrices,  respec- 
tively, and  x(t)  is  an  n-dimensional  state  vector,  u(t)  is  an  r-dimensional 
control  vector,  and  y(t)  is  an  m-dimensional  output  vector,  respectively. 
Further,  we  assume  that  the  system  (2.1)  is  completely  observable.   The 
optimal  output  regulator  problem  is  then  to  determine  the  control  u(t)  which 
minimizes  the  quadratic  cost  functional, 

v,  =  |yT  (O  Fy(tJ  4  f  (yT(t)  Qy(t)  +  uT(t )Ru(t ))at 

tQ  (2.2) 

where  P  and  Q  are  m  x  m  symmetric  positive   semidefinite  matrices   and  R  is 
an  r  x  r  symmetric  positive  definite  matrix.      Substituting  y(t)   =   Cx(t) 
into    (2.2),    V     is   rewritten  as 

Vl   =  2   xT(tf}    cTpCx(tf)    +  \     1      (xT(t)CTQC  x(t)    +  uT(t)Ru(t))dt 

"  t0  (2.3) 

T        T 
Since  the  system  (2.1)  is  completely  observable,  C  PC  and  C  QC  are  symmetric 

and  positive  semidefinite  [l]. 


The  optimal  feedback  control,  therefore,  is  given  by 

u*(t)  =  -R_1BTK(t)x(t)  -  (2.10 

where  K(t )  is  an  n  x  n  symmetric  positive  definite  matrix  which  is  the 


solution  of  the  matrix  Riccati  equation 

K(t)  +  K(t)  A  +  ATK(t)  -  K(t)BR_1BTK(t)  +  CTQC  =  0         (2.5) 
with  the  boundary  condition 

K(tf)  =  CTPC  (2.6) 

Furthermore,  the  optimal  trajectory  is  the  solution  of  the  system  of 
differential  equations 

x(t)  =  (A  -  BR~1BTK(t))  x(t)  (2.7) 

x(tQ)  =  xQ 

The  minimum  cost  is  given  by 

V1*  =  |xT(t)  K(t)  x(t)  (2.8) 

In  addition  to  the  assumption  of  complete  observability  we  further 
assume  that  the  system  is  completely  controllable.   Setting  P  =  0  and 
t   =  °°,  the  output  regulator  problem  turns  out  to  be  the  following: 
minimize 

V  =  |  /~(yT(t)Qy(t)  +  uT(t)Ru(t))  dt  (2.9) 

subject  to 

x(t)  =  Ax(t)  +  Bu(t) 

y(t)  =  Cx(t)  (2.10) 

x(0)  =  xQ 
Then  the  optimal  control  is  given  by 

u(t)  =  -R_1BTK  x(t)  (2.11) 

where  K  is  the  n  x  n  symmetric  positive  definite  matrix  which  is  the  solu- 
tion of  algebraic  Riccati  equation 

KA  +  ATK  -  KBR_1BTK  +  CTQC  =  0  (2.12) 

The  optimal  trajectory  is  the  solution  of  the  system 


i(t)  =  Gx(t) 


-1  T 
G  =  A  -  BR  B  K 


(2.13) 


:(0) 


^0 


Note  that  for  an  optimal  control  system  the  matrix  G  is  diagonable  and  has' 
eigenvalues  with  negative  real  parts.   Therefore,  the  symmetric  positive 
definite  solution  K  of  (2.12)  must  be  such'  that  G  has  eigenvalues  "with 
negative  real  parts . 

Kalman  [2]  has  shown  that  if  the  system  (2.1)  is  completely  observable 
and  completely  controllable,  then 

lim  K(t  )  =  K  (2.11+) 

t  -*» 

f 


We  shall  observe  this  equation  in  the  following, 
Let  K(t)  =  Z(t)X-1(t) 


(2.15) 


where  X(t)   and  Z(t)   are  n  x  n  matrices    and  X(t )    is   non-singular.      Substituting 
K(t)   =    (2(t)   -  K(t)   X(t))   X_1(t)  (2.16) 

into    (2.5)    and  setting 

X(t)   =  AX(t)   -   BR_1BT   Z(t)  (2.17) 

we  obtain 

Z(t)   =  -CTQCX(t)   -  ATZ(t)  (2.18) 

Therefore    (2.5)  becomes    a  system  of  2n-dimensional   linear  homogeneous 

differential   equations 


"x(t)  i     r 

_Z(t)  J       [. 


-1  T 
-BR  B 


T        T 
-C  QC    -A 


X(t) 
Z(t) 


(2.19) 


Since  we  are  interested  in  increasing  the  time  t,  let 


t  =  tf  -  t 


(2.20) 


Then 


X(t) 

Z(t) 


=  ¥ 


X(t) 
Z(t) 


(2.21 


where 


W  = 


-A 

T 
C  QC 


-1  T 
BR  B 


with 


X(0)  =  I 


Z(0)  =  CXQC  (2.22) 

Therefore,  the  solution  will  be  given  by  the  following  matrix  exponential 


X(t) 
Z(t) 


=  e 


Wx 


T 
C  QC 


(2.23) 


Consequently,  the  solution  of  the  matrix  Riccati  equation  will  be  obtained 
by  (2.15)  as  t->  °°.   In  the  next  section,  we  shall  show  that  two  matrices 
X(t)  and  Z(x)  are  obtained  from  the  similarity  transformation  which  reduces 
W  into  the  Jordan  canonical  form.   Throughout  the  remaining  discussions, 
the  increasing  time  variable  t  will  be  used  instead  of  x. 


3.   SOLUTION  OF  THE  MATRIX  RICCATI  EQUATION 


We  begin  with  two  lemmas  for  our  main  result. 

Lemma  1.   The  eigenvalues  of  the  matrix  W  of  (2.21)  are  symmetric 
with  respect  to  the  imaginary  axis. 


Proof.   Let  I  be  a  2n  x  2n  matrix  given  by 


xi  = 


-I 
0 


where   I   is   an  n  x  n   identity  matrix.      Then   it    is  obvious  that 

I  -1  =   I  T 
1  1 

T  T 

I  WI        =   -W 


(3.1) 


("3. 2) 
(3.3) 


Since  the   eigenvalues   of  a  matrix  are  unchanged  by   either  similarity  trans- 
formations  or  transposing,   the   eigenvalues   of  W  occur   in  pairs  with  opposite 
signs . 

Lemma  2.      Let 


T  = 


All       A12 


A  A 

LA21        A22J 


(3.M 


where  A      ,   A      ,  A      ,    and  A        are  n  x  n,  m  x  m,   n  x  m,    and  m  x  n  matrices, 
respectively.      Assume  that   A^      is  non-singular.      Then 


„-l 


T  T 

11  12 


T  T 

L     21  22J 


(3.5) 


exists  iff  A   -  A   A     A.^   is  non-singular,  and  is  given  by 


T__  =  A,,"1  -  -.  .   u  A  „  T. 


"11     11 


1   ^2  T22  A21  All  ' 


-1 


T    =  -A    ""  A    T 
12    All   A12  X22 


3.6) 


T21   "T22  A21  ^1 


-1 


T   =  (A   -  A   A 
22    VA22   A21   11 


-1 


*u.» 


-1 


where  submat rices  T^,  T^,    T±2,    and  T21  are  the  same  sizes  as  A 
A22'  A12'  and  A21'  resPectively- 
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-1 


-1 


Proof.   Suppose  Ag2  -  Agl  A^    AL2  is  non-singular.   From  TT   =  I, 


All  Tll  +  A12  T21  '  :„ 


All  T12  +  ^2  T22  =  ° 


A21  Tll  +  A22  T21  "  ° 


(3.7) 


A21  T12  +  A22  T22  =  Im 


-1 


Premultiplying  the   second  equation  by  A       A         "   and  subtracting  from  the 


fourth, 


T22    (A22  "  A21  All  '"  A12}  ' 


And  hence 


-1 


T12    ^1   h2   T22 


Similarly,  from  the  first  and  third, 
T21  =  "T22  A21  ^Ll 


Tu  ■  An_1  -  All  *  v 


"I        m  "I 

^2  T22  A21  All 


Hence  T   is  a  right  inverse  of  T.   In  a  similar  way  we  can  show  that  T 


is  also  a  left  inverse  of  T.   Conversely,  suppose  T  is  non-singular.   Then 

-1      . 


0  ?     det 


^.1     A12 


A  A 

A21     A22J 


=  det 


\i  0l 


dl       m 


]  F : 


12 


°        A22"A21  \l 


-1 


12  J 


-1 


"  det    <All'    det(A22  "  A21  *U  "  A12) 

Therefore,   A       -  A       A  A        is  nonsingular. 


(3.8) 


Assuming  that  W  has  distinct  eigenvalues,  we  now  turn  to  constructing 
the  similarity  transformation  with  which  ¥  of  (2.2l)  is  reduced  to  the 


diagonal  form.   By  lemma  1,  W  has  n  eigenvalues  with  positive  real  parts 
and  n  eigenvalues  with  negative  real  parts.   If  we  construct  the  similarity 
transformation  S,  whose  columns  are  the  right  eigenvectors  of  W,  such  that 
the  first  n  columns  are  the  right  eigenvectors  of  W  corresponding  to  the 
n  eigenvalues  of  W  with  positive  real  parts,  then 

rx„ 


S_1¥S  = 


X. 


(3.9) 


2n 

where,   of  course,   the  rows   of  S        are  the  left   eigenvectors   of  W. 
O'donnell    [5]   and  Potter    [6]   have   shown  that   if  we  partition  S   into   four 
n  x  n  submatrices, 


S  = 


Sll        S12 


.S21        S22J 


(3.10) 


then  the  solution  of  the  matrix  Riccati  equation  is  given  by 

K  -  S21  hl^ 
provided  that  S   is  non-singular. 


(3.11) 


We   shall  now  remove  the   restriction  that  W  has    distinct   eigenvalues   and 
present   an  alternative  proof  to  that   of  Martensson    [7]   to   show  that    (3.1l) 
holds   for  the   case   in  which  W  has   nondiagonal  Jordan  form. 

Let  us   consider  the   general  Jordan  form  J  with  the   similarity  transfor- 
mation S,   which  will  "be  constructed  later, 
S_1¥S  =  J 


where,      J  = 


*! 


(3.12) 


m, 


and  in  which 


m. 


A.    1 

1 


A.    1 
l 


\  A. 

l 


,  (m.  x  m.  ) 


(3.13) 


If  we  denote  S  by  its  column  vectors  s, .  ....  s^  ,  then 

1        2n 

¥  Lsl5  s2,  ...,  s2n]  =  [s1,  s2,  ...,  s2n]  J 


(3.1*0 


For  any  Jordan  block  J      ,    if  m.    >  1  then,   dropping  subscripts  of  s   for 

i 
convenience. 


Ws(l)=s(1)A 


ws  =   s  +   A  .    s 

1 


(3.15) 


or 


Ws 


(m. )   .       (m.-l)  (m.  ) 

1=S1  +A.S1 


(¥  -    A.I)    s(l)    =   0 

l 

(¥  -    A.I)2   s(2)    =   0 


(3.16) 


/  sm.       (m. ) 

(W-A.Ijis      i     =0 


Thus,  s    ,  1  <  j  <  m. ,  is  a  principal  vector  of  degree  j  associated  with 

(l  )  (m    ) 

A..   Since  S  is  nonsingular,  the  m.  principal  vectors  s   ,  ...,  s  i 

are  independent.   If  we  consider  all  the  Jordan  blocks  J   ,  1  <  i  <  £,  we 

m.    ~ 

l 

have  the  2n  principal  vectors  s  ,  s  ,  ...,  s   which  constitute  the 
similarity  transformation  S  in  (3.12).   Further,  we  rearrange  the  columns 
of  S  such  that  the  first  n  vectors  are  the  principal  vectors  corresponding 
to  the  n  eigenvalues  with  positive  real  parts  of  W. 
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Let  us  now  split  J  of  (3.12)  intc 


J  = 


(3.17) 


where  J     and  J     are  diagonal  with  the  eigenvalues   of  W  with  positive  and 

negative  real  parts  respectively,    and  E     and  Ep   are  matrices  with  one's  on 

the  super-diagonals  and  zeros   elsewhere.      To   associate  with    (3.6),   let   us 

partition  S  and  S  as 


s 

S  „ 

11 

12 

s  = 

s 

s„„ 

21 

22 

^ 

T 

T 

0-l 

11 

12 

s  = 

T 

T 

21 

22 

(3.18) 


where  all   submatrices   are  n  x  n.      Note  that,  by  lemma  2,    it   is  necessary  to 
assume  that   S        and  S22  -   S        S     ~     S        are  nonsingular. 


Then 


Wt  SJS-1t 

e        =  e 


-  Se^s"1 


sll 

S12 

f 

S21 

S22 

Jit    n 

e  1        0 


0        e  <= 


At 


E2t 


T  T 

11  12 


T  T 

21  22 


S11eJlteEltT11  +   Sl2eWS\       ^eWl^  +   ^e^ 


J.t  E,t 


J«t  E„t 


J.  t  Entr 


S-_e  1  e  1  T       +    S./2V2X,        S_eulWT       +    R„ed2V2lT 


21 


J^t  E„tr 


"11  22 


21  21 


12  22 


"22 
(3.19) 
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Since  the  elements   of  J     have  negative  real  parts,    e  2     -*•  0  as  t  ->  °°     and 
thus   the  terms  having  e   2     in    (3.19)   vanish   as  t  -*  °°   .      Therefore, 


X(t) 
Z(t) 


=  e 


Wt 


T 
C  QC 


ri 


T 
C   QC 


S1:L   eV   eV    (T12   +  T^   CTQC)' 
S21   e*7!*   eV    (T       +  T        CTQC ) 


(3.20) 


Therefore  the   solution  of  the  matrix  Riccati   equation  is   given  by 
K  =  lim     Z(t)   X_1(t) 

=  lim   [S        eV   eV    (T-.+T       CTQC)  ]  [S     eVeV  (T     +T     CTQC)  I"1 


■  S21   Sll 


-1 


(3-21) 


provided  that  the  indicated  inverse  exists.  However,  we  shall  show  that 
only  the  nonsingularity  of  S  and  S  -  S  S  S  ?  is  needed  for  the 
stationary  solution  of  the  matrix  Riccati  equation.   From  (3.12), 


T    T 
11   12 


T    T 
21   22 


-1  T 
-A   BR  B 


T       T 
C  QC    A 


Sll   S12 


S21   S22 


=  J 


(3.22) 


Equating  the  (2,1 )  elements  on  each  side  of  (3.22), 


(-T21  A  +  T22  CTQC)  S1X  +  (T21  BR_1BT  +  T^  AT)  S21  =  0 


12 


or 


"T21  AS11  +  T22  ^  S21  +  T21  BR_1bT  S21+  T22  ^^  Sll  =  °     (3'23) 
Substituting  T   =  -  T22  S21  S^"  into  (3.23), 

T22S21Sll"1+  T22ATS21  "  T22S21SirlBR"lBTs21  +  T22cT«°  Sll=°  (3  ^ 

Premultiplying  by  T     and  postmulti plying  by  S    , 

S21S11_1a  +  aTs21S11_1  "  S21S1i"1bR"1bTS21S11_1  +  ^  =  °    (3'25) 
which  is  to  be  proved. 
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k.      COMPUTATIONAL  METHODS  AND  NUMERICAL  RESULTS 

The  major  three  steps  for  computing  the  solution  of  the  matrix  Riccati 
equation  are  as  follows: 

1)  Compute  eigenvalues  and  eigenvectors  of  W, 

2)  Compute  principal  vectors  of  W  if  W  is  defective, 

3)  Compute  S   S  "   from  S. 

For  computing  the  eigensystem,  we  used  two  methods:   (i)  the  QR 
algorithm  [9]  with  an  inverse  iteration  routine  [10],  and  (ii)  Jacobi-like 
method  [ll]  for  finding  the  eigenvalues  and  eigenvectors.   Their  perfor- 
mances are  satisfactory  and  two  numerical  test  examples  are  given  in  Test 
Result  1  and  Test  Result  3.   When  the  matrix  W  is  nondiagonable,  Test  Result 
2,  the  first  method  requires  the  evaluation  of  the  principal  vectors  as 
well  [12].   This  poses  some  numerical  difficulties  as  we  are  essentially 
trying  to  solve  the  ill-conditioned  system  of  equation 

(W  -  Al)k  x    =  x.,  k  >  1 
where  ^  is  an  approximation  of  true  eigenvalues  A.   Thus,  the  Jacobi-like 
method  is  more  attractive  in  this  case  for  providing  an  approximate  eigen- 
system that  is  suitable  for  engineering  purposes. 

Appendix  I  contains  the  listings  of  the  ALGOL  programs  used  on  the  B-6500 
Appendix  II  contains  the  listings  of  two  ILLIAC  IV  programs,  namely,  the 
inverse  iteration  routine  and  a  routine  for  solving  a  system  of  complex 
linear  equations,  both  written  in  GLYPNIR.   ILLIAC  IV  programs  for  the  QR 
algorithm  [13,  lM  and  the  Jacobi-like  algorithms  [15,  l6]  are  already 
available. 
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Test   Result   1, 


A  = 


1 

0 

0 

0 

0 

1 

0 

-1 

0 

0 

0 

0 

-1 

0 

0 

0 

0 

1 

0 

-1 

0 

0 

0 

0 

-1 

p  = 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

B  = 


1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

Q  = 


0 

0 

0 

0 

0' 

0 

10 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

10 

0 

.0 

0 

0 

0 

0 

c  = 


1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

R  = 


10  0 
0  10 
0        0        1 


The  numerical   example  was  tested  "by  the   QR  algorithm  and  the  inverse 
iteration  method.      The  ten  distinct   eigenvalues   of  W  are 

±   1, 

±  1.7287601*77185  ±  1.577533767573i, 

±  1.35319578U0U6  ±  l.l537^989928Ui. 
We  choose  the  five  eigenvectors  corresponding  to  the  five  eigenvalues 
with  positive  real  parts  and  computing  K  =  S   S  ~   from  them  gives  the 
solution  K  as 
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1.262782609 
2.  i+9^009759 

-0.819173651 

0.668267901 

-O.UU3608958 


2.1+9^009759 
7.1+351+51161+ 

-1.8257^1858 
1.122910U32 

-0.668267901 


-0.819173651 
-I.8257U1858 
1.6383U7303 
1.8257^1858 
-0.819173651 


0.668267901 
1.1229101+32 
1.8257^1858 
7. 1+351+ 511 6k 
-2.1+91+009759 


-0.1+1+ 3608958 
-0.668267901 
-O.819173651 
-2.1+91+009759 
I.262782609 


All  the   eigenvalues   of  G  =  A  -  BR     B  K  have  negative   real  parts   and  hence 
the   above   solution  K  is   acceptable. 


Test   Result   2. 


A  = 


■3        2 
■2        1 


B  = 


P=Q= 


0        0 
0        0 


The  mat 

rix 

3 

-2 

0 

0 

2 

-1 

0 

1 

W  = 

0 

0 

-3 

-2 

0 

0 

2 

1 

c  = 


R  = 


1        0 
0        1 

CO 


is   defective,   that    is,   W  has    four  eigenvalues   ±1    (double)   and  each  +1   and 
-1   correspond  to  a  quadratic   elementary   divisor,    respectively.      This 
example  was  tested  both  by  Jacobi-like  method  and  by  a  combination  of  the 
QR     algorithm  and  the   inverse   iteration  method. 
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The  Jacobi-like  method  gives   K  as 
-10.000008*+  6.000001+3 

6.0000031+  -k.  000002T 

and  the  QR  and  inverse  iteration  method  give  K  as 
-9. 9999917  5.99999^7 

5.9999953         -3.9999983 
The  exact   solution  is 

"Vlo      6 
6     -1+ 

For  computing  the  principal  vector  "by  the  QR  and  inverse   iteration,  we 
chose  the  perturbation   for  the  true  eigenvalues   as  much  as  10 


Test   Result   3- 


A  = 


6 

h 

1+ 

1 

1+ 

6 

1 

k 

1+ 

1 

6 

k 

1 

1+ 

k 

6 

B  = 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

C=R= 


1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

P=Q= 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

The  exact  eigenvalues  of  W  are 

±15,  ±5,  ±5,  ±1 
Since  the  matrix  W  is  derogatory,  +5  and  -5  correspond  to  a  two-dimensional 
subspace,  respectively.   The  computed  eigenvalues  by  Jacobi-like  method 


are 
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Ik. 99999999930 

^.9999999991+9 
-h. 99999999982 

0.99999999997 


-1^.9999999981+ 
U. 99999999981 
j+.  999999999)48 
-0.99999999995 


The   computed  solution  K  is   given  "by 


■1.999999999811 

1.99999999988U 

2.000000000276 

-2.00000000036U 


1.9999999999U2  1.999999999913 
-2.000000000015  -2.000000000000 

-2.0000000001+07  -2.000000000378 
2.000000000509   2.0000000001*80 


-2.000000000015 
2.000000000087 
2.000000000U80 

-2.000000000568 


in  which  the  absolute  error  is  less  than  10 


-9 
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APPENDIX  I 


ALGOL  PROGRAMS 
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t  TWO  ALGot  PROCEDURES'  RlCCATU  AnO  RICCATI2.  FOR  COMMUTING  THE  OOOOOlOO? 

f  SOLUTION  OF  THE  MATRIX  RICCATI  EQUATION  ARE  LISTED  BtLOW*  FOR  00000200? 

(  «RICCATlU.  THE  RROCtOURE  "EIGEN*  (BY  EBERLEIN  AND  BUQThNoyD)  WAS  00000300? 

I  USEU  WHILE  A  COMBINATION  OF  "HQR"  (BY  MARTIN.  ET  AL)  AND  00000400? 

I  MNVITErATION"  WAS  APPLIED  TO  "RICCATI2* f  LINEAR  EqUATIUn  SOLVER  00000500? 
%    (BY  BOWDlER.  ET  AL)  WAS  ALSO  USED  TO  PROVIDE  THE  INytRSE  OF  A  MATRIX   00000600? 

«  ANO  TO  CALCULATE  THE  PRINCIPAL  VECTORS  FROH  ( ( W-LAMBUAM  >**K>*X«0.  00000700? 

00000800? 

00000900? 

PRDCEUURL  RlCCATIl(N.L»M»A»B#C»P»O.R»THX»W»K)l  00001000? 

VALUE  N.L.M.TMXI  00001100? 

INTEGtR  N,L.M.TMX|  00001200? 

ARRAY  A.B.C.P»a«R»W»K[l»ll|  00001300? 
COMMENT  SOLVES  ThE  MATRIX  RlccATI  EQUATION                              .  OOOOUOO? 

KA*(AT)K-KR(RINV)(BT)K*(CT)QC«o,  00001500? 

A  Is  N*Nt  R  IS  N*L»  C  IS  M*N»  P  AND  0  ARE  M*M.  R  IS  L*L  OOOOUOO? 

MATRICES*  RESPECTIVELY*  FORM  2N*2n  MATRIx  W  ANU  COMPUTE  00001700? 

EIGENVALUES  ANU  EIGENVECTORS  OF  W*  WHERE  FOUR  N*N  BLOCKS  OOOOUOO? 

Op  w  ARE  (1»1>B*A»  (1#2)«B(RINV)<BT)»  < 2» 1 )■( C I )«C»  ANO  00001900? 

<2»2>"<AT>.  REARRANGE  EI GEN VEcTOK C Pr INC  I PAL  vEcTOR>  HATRlX  00002000? 

SO  THAT  ThE  rlHST  N  COLUMNS  CORRESPOND  TO  TO  THE  EIGENVALUES  00002100? 

WlTH  P0SITIVE  REAL  pARTS.  dENoTINq  ThE  UPP^R  HaLf  Op  N  COLUMN  00002200? 

VECTORS  BY  SMTXl  ANO  THE  LOWER  HALF  BY  SMTX2.  I  HE  SOLUTION  00002300? 

K  IS  THEN  GlvEN  BY  K»C SMTx2) ( SMTxl INv ) .  THE  PRUCEDURE  00002a00? 

"plGEN"  IS  USED  ^OR  EIGEnPrUBLEM  0F  Wi  0000250O? 

00002600? 

PEGIN  COMMENT  plRST  DECLARE  PROCEOuRESl  00002700? 

00002800? 

PROCEUURE  SETUP(N.L»M»A»B»C»Q»RINV.W)|  00002900? 

VALUE  N.L.M,   INTEGER  N»L»M>  00003000? 

REAL  ARRAY  A • B# C » Q . R I N V # W (  1  »  1]  I  00003100? 

COMMENT   THE  PROCEDURE  SETS  UP  ThE  2N*2N  MATRIX  W.  rInV  IS  INVERSE  Op  RI00003200? 

BEGIN  OOOO33OO? 

REAL  ARRAY  R I NVBTU  I  L»  1  I  N] .  QC  [  1  I  M.  1  I  N  ]  |  00003*00? 

REAL  SUM)  00003500? 

INTtGER  I.J.KI  00003600? 

FOR  It"i  STEP  1  UNTIL  N  DO  00003700? 

TOR  Jl«l  STEP  1  UNTIL  N  DO  00003800? 

8LGIN  00003Q00? 

W(I.J1I.-A(I,J]|  00004000? 

W(N*I,N*J]l«AtJ.I]»  00004100? 

ENUl  00004200? 

FOR  Il«l  STEP  1  UNTIL  L  00  00004300? 

FOR  J»«l  STEP  1  UNTIL  N  DO  00004*00? 

BEGIN  StlHl.O.Ol  00004500? 

FOR  K i-l  STEP  1  UNTIL  L  DO  00004600? 

SuM|«SUM*RINV(I»K]*B[J»KJI  00004700? 

RINvBTd.JJl-SUMl  00004800? 

ENOl  00004900? 

TOR  I|»l  STEP  1  UNTIL  N  00  00?Sc°207 

FOR  Jl-1  STEP  1  UNTIL  N  00  00005100? 

BttilN  SyMl.O.Ol  00005200? 

FOR  Kill  STEP  1  UNTIL  L  00  22SS!3221 

SuH|«SUM*B[I,K]*RlNVBTtK»J]|  0000*400? 

"CI.N*J]l«SUMl  00003500? 


END| 


00005600? 


FOR  If .1  STEP  1  UNTIL  M  DO  00005700? 
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FOR  Jial  STEP  1  UNTIL  N  00  000058001 

BtGlN  SljMlaO.Oj  000059001 

FOR  Kt«l  STEP  1  UNTIL  M  00  000060001 

SuMI«SUM*QtI»KJ*C[K»J]l  000061001 

OCC I • J] laSUMl  0000*2001 

EN0|  000063001 

FOR  I|«l  STEP  I  UNTIL  N  DO  00006*001 

FOR  J  is  1  STEP  i  UNTIL  N  00  000065001 

BtGlN  $(jMl«0,Ot  000066001 

FOR  Kf«l  STEP  1  UNTIL  N  00  000067001 

SijMI"SUM*CCKtlJ*ac[K»JJI  000068001 

W[N*I,J)|«SUMI  000069001 

ENOi  000070001 

ENO  OF  SETUPj  000071001 

000072001 

COMMENT  INSERT  THE  PROCEDURES  INNPROD. DECOMPOSE.  AND  ACCSOLvEi           000073001 

oooofaoo? 

PROCEUURE  INVERSE(N»A.X)|  000075001 

VALUE  N|  INTEGER  N,  000076001 

DOUBLE  ARRAY  A.XU.ljf  0000/7001 

COMMENT  TmE  PROCEDURE  COMPUTES  THE  INVERSE  OF  AN  N*N  RtAL                000078001 

UNSyMMrTRIc  MATRIX  A,  USINq  The  PROCeOUReS  "DrC0MPU$E«  AND  0000?900l 

"ABSOLVE",  WRITING  ON  X|  000080001 

REGIN  000081001 

INTtGER  UJ.IT.D2l  000082001 

REAL  Oil  000083001 

ARRAY  AA.RB.BCllN.l |N]*PC1 INll  00008*001 

FUR  Il«l  STEP  1  UNTIL  N  DO  000085001 

FuR  Jl«l  STEP  1  UNTIL  N  Do  000086001 

IF  J«I  THEN  BtI.J]»«1.0  ELSE  BCI»J]l«0.0>  000087001 

FUR  I,«i  STEP  1  UNTIL  N  DO  000088001 

FUR  J i«l  STEP  1  UNTIL  N  DO  000089001 

AAC I.J] laA[I,J]|  000090001 

DtC0Mp0sE(N»AA,0l.D2.P)l  000091001 

AtCSULVE(N.N.A.AA.P.B.X.BB»IT)|  000092001 

ENO  INVERsEl  000093001 

00009*001 

comment   Insert  the  procedure  eigeni  000095001 

000096001 

PROCEUURE  MuLT(Nl.N2,NJ»AiB.C)l  000097001 

VALUE  N1»n2.N3|  iNTEfiEH  Ni»N2.N3l  000098001 

DOUBLt  ARRAY  A.B.Ctl.UI  000099001 

REGIN  000100001 

INTtGER  I.J.KI  000101001 

REAL  SUM1  00010200? 

FOR  It«l  STEP  1  UNTIL  Ni  DO  00010300? 

FOR  J i ■  1  STEP  1  UNTIL  N3  00  00010*00? 

BEGIN  SuH|aO.Ol  00010500? 

FUR  k,"1  STEP  1  UNTIL  N2  00  00010600? 

SUMI»SUM*A[I»K)*BCK,JJ|  00010  700? 

C I  I . J ] I.SUMI  00010800? 

END*  00010900? 

END  MULTl  00011000? 

000111007 

ARRAY  HiNV»SMTXl,SMTX2UlN.llN].wW.TEMP,VECT(llN*N.l!N*NJ|  00011200? 

INTLGER  ArRAy  SEARCHlHNll  00011300? 

iNTtGEH  N2»I.J.InDEX»  00011*00? 
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N2I*N*N|  OoOlHOO? 

INVtRS£(L.R.RlNV)l  OOOilloo? 

SETUP<N,L,M»A«B»C»Q»RINV»W)J  0001  If 00? 

FOR  I • » 1  STEP  1  UNTIL  N2  DO  00011800? 

FOR  Jl«l  STEP  1  UNTIL  N2  00  00011900? 

TtMRCi.jji.wti.jji  00012000? 

EIG6N(N2»TMX.TEMP.VECT)I  00012100? 

COMMENT  SEARCM  N  COLUMNS  OF  VECT  CORRESPONDING  TO  N  6lGENVALUES  00012200? 

WITH  P0SITIV6  REAL  PARTS.  SNTX1  ANO  SMTX2  ARt  IT*  UPPER  00012300? 

AND  LOWER  HALVES*  RESPECTIVELYi  00012400? 

IN0LX|«01  00012500? 

FOR  Il«l  STEp  1  UNTIL  N2  DO  00012*00? 

If  TEmP[I»I]  GTR  0  THEN  00012700? 

3LGIN  00012800? 

IN0EX|.INDEX*1  ,  00012900? 

SEARCHtlNDEXjl.il  00013000? 

ENDi  00013100? 

FOR  Ji.l  STEP  1  UNTIL  N  00  00013200? 

FOR  It.l  STEP  1  UNTIL  N  00  00013300? 

REG1N  00013400? 

SMTxlt Ii J) i.VEcTtl'SEARCMC Jll  |  00013500? 

SMTX2tI.J]l"VEcTtN*I»SEARcHtJm  00013600? 

END!  00013^00? 

INVtRSE(N.SMTXl.TEMP)l  00013800? 

MULICn»N»N»SMTx2»TEMP»K)|  00013900? 

End  riccatH»  ooomooo? 

ooohioo? 

00014200? 

PR0CEUURE  RiCCATi2(N.|.»M»A#BtC»P»8,R»W»K.DEfECT»OERoGAlE)l  0001 4 300? 

VALUE  N.L.M|  00014400? 

INTEGLR  N.L,M|  00014500? 

ARRAY  A.3.C.Pti]fRiMfK{  lilJl  00014600? 

INTEGtH  ARRAY  DEFECT .D6HUGATEI 1  J  I  00014700? 

COMMENT  SnLvES  THE  MATRIX  RICCATI  EQUATION  00014800? 

KA*(AT)K-KR(KINV)(BT)K*(CT)aCO.  00014900? 

A  IS  N*N.  h  IS  N*L»  C  IS  M#N.  P  AND  Q  ARE  M*M,  R  IS  L*L  00015000? 

matrices.  resptctively.  form  2n*2n  matrix  w  anu  compute  00015100? 

eigenvalues  anu  eigenvectors  of  w,  where  four  n*n  slocks  00015200? 

or  w  are  c i • 1 >"-a.  (1.2).b(rinv)(bt)»  ( 2» 1 >■< c  i  >0c»  and  00015300? 

C?»2)»(AT).  REARRANGE  E1GEnVECT0R(PrINCIPAL  VECTOR)  MATRIX  00015400? 

SO  THAT  THE  rlKST  N  COLUMNS  CORRESPOND  TO  TO  TM6  EIGENVALUES  00015500? 
«ITh  POSITIVE  HEAL  PARTS.  DENOTING  THE  UPPER  H*LF  "F  N  COLUMN    00015600? 

VECTORS  BY  SHTXi  AND  THE  LOWER  HALF  BY  SMTX2»  THE  SOLUTION  00015700? 

K  Is  THEN  GIVEN  BY  K.( SMTX2) ( SMTX1 INV ) .  THE  PRUCEDURES  00015800? 

"mShLD"'  "hOr"»  aN0  "INVITEHATlON"  ARE  USED  TOM  EI<»ENPRORlEM  00015900? 

Op  M|  00016000? 

00016100? 

BEGIN  COMMENT  FIRST  DECLARE  PROCEDURES!  00016200? 

00016300? 

PROCEDURE  SE!TUP(N.L»M.A»B#C»a»RINV.W)|  OOOI64OO? 

VALUE  N.L.Mi   INTEGER  N»L.M»  00016500? 

RF.AL  ARRAY  A » B » C » Q . R I N V • W t  1  »  1  J I  00016600? 

COMMENT  THE  PROCEDURE  SETS  UP  ThE  2N*2N  MATRIX  W.  RINV  IS  INVERSE  Op  Rl  00016700? 

PEGIN  00016800? 

REAL  ARRAY  R I N VBT t 1  I  Li  1  I N J . QC t 1  I H . 1  I N J  I  00016900? 

REAL  SUmi  00017000? 

iNTtGER  I.JfKI  00017100? 
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FOR  II" 

for  Ji« 

BtGIN 

W[I 

W[N 

ENOI 

FoR  Il« 

FOR  Ji« 

BtGIN 

FOR 

S 

RIN 

ENO| 

FOR  I|» 

FOR  Jl" 

BLGIN 

FOR 

S 

wc  I 

ENOi 
FOR  Iia 
FOR  J|" 

BtQlN 

FOR 
S 

act 
enOi 

FOR  Il» 
FOR  J|« 

BLGIN 

FOR 

S 

NCN 

END| 

END  OF  SE 


1  STCP  1  UNT 
1  STEP  1  UNT 

, J] I«-At I ,J1 
♦I.N*J]I«A[J 

1  STEP  1  UNT 
1  STEP  1  UNT 
SUM|«0»0| 
K|»l  STEp  1 

(jM|«SUM*RINv 
VBTC I* J] l«SU 

1  STEP  1  UNT 
1  STEp  1  UNT 

SuMl«0,0| 
K|«l  STEP  1 

UM,«SUM*BtI. 

.N  +  J]  |"SUm» 

1  STEP  1  UNT 
1  STEP  1  UNT 
Sum«"0«0| 

K|»l  STEP  1 
UM|»SUM*Q(  I, 

I# jl laSUMl 

1  STEP  1  UNT 
1  STEP  1  UNT 
SUMl«0,Ol 
K |«1  STEP  I 

UM|«SUM*CtK. 
♦I, J] l"SUMl 

TUPI 


IL  N  DO 
IL  N  DO 

I 
•  I  Jl 

I*.  I    Do 
II  N  DO 

UNTIL  L  DO 
f  lfK]*B[Jf K]l 
Ml 

IL  N  DO 
IL  N  00 

UNTIL  L  DO 
Kl*RINVBTCKf JJ| 


IL  M  00 

IL  N  DO 

UNTIL  M  DO 
K]*C[Kf J] I 


IL  N  DO 
IL  N  DO 

UNTIL  M  DO 
I]*UC[K»JJI 


COMMENT  INSERT  THE  PROCEDURES  INNpROD. DECOMPOSE.  AND  AtCSOLvEl 

PROCEDURE  INVERSE(N.A.X)! 

VALUE  N|  INtEqEH  N| 

DOUBLE  ARRAY  A.XU.lJI 

COMMENT  ThE  PROCEDURE  COMPUTES  THE  INVERSE  OF  AN  N*N  RtA|_ 

unsymmetric  matrix  a.  using  the  procedures  "DECQmpUse" 

"ACCSOLvE".  writing  ON  Xl 
REGIN 

INTLGER  I.J. IT. 021 
REAL  ol| 

ARRAY  AA.BB.BlliN.l tNJ.Ptl|Nli 
FUR  Ii-i  STEP  1  UNTIL  N  DO 

FuR  J i "i  Step  i  until  n  d0 

IF  J«I  THEN  B[I.J]I»1.0  else 
FUR  I,«i  STEP  1  UNTIL  N  DO 
FUR  Ji-i  STEP  1  UNTIL  N  DO 

AAlIf J] l«A[I, J]| 

DECOMPOSE (Nf AA.0l.D2. P)l 

ACCSOLVE(N.N.A.AA.P.B.X.BB.IT)i 
END  INVERsEi 


AND 


B[  If  J]  i»0.0J 


00017200? 

000173001 

00017400? 

00017500? 

00017600? 

00017700? 

000179007 

000179007B 

00018000?) 

00018100? 

00018200? 

00018300? 

00018400? 

00018500? 

00018600? 

00018700? 

00018800? 

00018900? 

00019000? 

000i9i00? 

00019200? 

00019300? 

00019400? 

00019^00? 

00019600? 

00019700? 

00019800? 

00019900? 

00020000? 

00020100? 

00020200? 

00020300? 

00020400? 

00020500? 

00020600? 

00020700? 

00020800? 

00020900? 

00021000? 

00021100? 

00021200? 

00021300? 

00021400? 

00021500? 

00021600? 

00021700? 

00021800? 

00021900? 

00022000? 

00022100? 

00022200? 

00022300? 

00022400? 

00022500? 

00022600? 

00022700? 

00022800? 
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comment  InSeRt  the  procedures  hshld  and  hori 

PROCEDURE  HuLT(Nl*N2.Ni*AtBtC)l 
VALUE  N1»n2,N3i  INTEQEK  Nl.N2.N3f 
DOUBLE  ARRAy  A»B*C[1.1)| 
BEGIN 

INTtOER  I.J»KI 

REAL  SUM! 

FOR  Ii«l  STEP  1  UNTIL  Nl  DO 

FOR  Jl«l  STEp  1  UNTIL  N3  DO 

BEGIN  SllM|«0,OI 

FUR  Ki«i  STEP  1  UNTIL  N2  DO 

SUM|«SUM*A[  I.K)*B[K. J] | 
C  t I . J  1  I -SUM  I 
ENOI 
ENO  MULTl 


PROCEUURE  lNVlTERATION(N»w»X»KK»P)| 
VALUE  N,  TNTEGER  N.KK| 

w.xci.n, 


'i 


DOUBLE  ARRAy 

array  pchi 
c0mment  s0lveS  wx 

u      MXTrI*.  wr 

THE  ElQENy 

N  Is  THE  0 

NEEoEO  for 

BEGIN 

INTEGER  I»J»D2»R 
DOUBLE  OI.MAXVAL 
DOUBLE  ARRAY  MH( 
LABEL  iNVlTiOUTl 
FOR  Ii"i  STEP  1 
FOR  J l -  1  STEP  I 

WW[ i, j]|«rf[  I.J 
DECUMpOsE(N.WW.D 
FOR  Il«i  STEP  1 

81 1.1 ] i.l .01 
KK I "1  J   R | ■  1  | 
INVITI 

ACCSOLVf;(N.R.  «»  W 
MAXVAl.l-0  j 
FOR    Ii"l     STEP    1 
If    ABS(xCI.l]) 
BEGIN 

MAXVALI"AB 
MAXjNOEXl* 
ENDi 

MAXVA|_|.X[MAXI 

FOR  I i » 1  STEP  1 

XI  I.  1  ]  I.X[  I.l  ] 
FOR  It»i  ST^P  1 
l>  ABS(BCI.I)- 
BtGIN 

if  kk  gtr  io 

FOR  J|«l  STE 
BtJ.l]fX[ 


■0  BY  INVERSE  ITERATION*  WHERE  *    IS  AN  N*N  REAL 
ITING  SOLUTION  ON  X,  THUS.  THE  PROCEUURE  COMPUTES 
ECTOHS  CORRESPONDING  TO  THE  GIVEN  EltiENv*LUES. 
RDER  UF  N  AND  KK  IS  THE  NUMBER  OF  iTtRATiONS 
ITERATING  EACH  ElGENvEcTORl 

.LiMAXlNOEXi 

I  |N.l|NJ.B.BBll|N.l  |U| 

UNTIL  N  00 
UNTIL  N  DO 

]l 

l.D2»P)» 
UNTIL  N  DO 


W.P.B»X.BB.L)J 

UNTIL  N  DO 
GTR  MAXVAL  THEN 

s(xri*n)i 
ii 

N0EX»1]I 

UNTIL  N  DO 

/MAXVALI 

UNTIL  N  00 

X[I.i]>  GTR  p-12  THEN 

then  go  to  outi 
p  1  until  n  du 

JM  1» 


00022900? 
00023000? 
00023100? 
00023200? 
00023300? 
00023400? 
00023500? 
00023600? 
00023700? 
00023800? 
00023900? 
00024000? 
00024100? 
00024200? 
00024300? 
00024400? 
00024500? 
00024600? 
00024700? 
00024800? 
00024900? 
00025000? 
00025100? 
00025200? 
00025300? 
00025400? 
00025500? 
00025600? 
00025700? 
00025600? 
00025900? 
00026000? 
00026100? 
00026200? 
00026300? 
00026400? 
00026500? 
00026600? 
00026700? 
00026800? 
00026900? 
00027000? 
00027100? 
00027200? 
00027300? 
00027400? 
00027500? 
00027*00? 
00027700? 
00027800? 
00027900? 

00028000? 
00028100? 

00028200? 
00028300? 
00028400? 
00028500? 


26 


KK  I  «KK*  1 1 

GO  TO  INVITI 
ENU| 

OuTl 

END  INVREsElTi 


INTLGER 
ARRAY  Te 

EI 
REAL  PER 
INTtGER 
N2l"N*N| 
INVtRSE( 
SETUP<N. 

FOR  I|«l 
FOR  Ji "1 

TLMRtl 
HSHLO(i  , 
HQR(N2»T 
FOR  I i - 1 
FOR  J  I ■ 1 

OKGWlI 
COMMENT 

FOR  JJ I ■ 

BEGIN 
FUR  It 
FUR  Ji 

FUR  It 

w[Ii 

I*    JJ« 
8EGI 

Fo 
My 

FO 

END  | 
IF  JJ. 

FOR 

«t 

INVITE 

FUR  It 

FUR  J| 

SMTX 

END  JJ» 
COMMENT 


I.J»JJ»CNT»n2»INDEX#ITC0UNT| 

Mp.0RGWtSMTX»VECT0Rti|N*Nil|N*N]»RlNVtSMTXlpSMTX2[l|N»ltN]. 

GR»EIGI»INTCH11«N*N]I 

TURB.MPLiERl 

ARRAY  SEARCHlllNll 

L»R#RINV)I 
L,M,A,B,C.Q.HINV,H)J 

STEP  1  UNTIL  N2  00 

STEP  1  UNTIL  N2  00 
• JllaMCl* Jll 
N-,TEMP)| 
EhP.EIGR.EIGI»CNT)| 

STEP  1  UNTIL  N2  00 

STEp  1  UNTIL  N2  00 
. J] laWCI, J]l 
COMPUTE  EIGENVECTORS 
1  STEP  1  UNTIL  N2  DO 


OR  PRINCIPAL  VECTORSl 


INDEX, «o 
FOR  U"l 
!►  ElG 
BLGIN 

INDE 
SEAR 
ENOi 


■1  STEP  1  UNTIL  N2  00 

■1  STEP  1  UNTIL  N2  DO 

J]  iBORGWt  I.JH 

■1  STEP  1  UNTIL  N2  DO 

I]«"WtI»i]-ElGR[ JJ)*PERTURB> 

DEFECT[JJ]  THEN 

N 

R  1 1*1  STEP  1  UNTIL  N2  00 

H[I»I]i«W(i,n*pERTURB*MpLIERI 

LT(N2»N2,N2»W»W,TEMP)I 

R  Ii.l  STEP  I  UNTIL  N2  DO 
R  Jlal  STEP  1  UNTIL  N2  DO 
*l I»J]««TEMP[  I»J)I 

DEROGATEtJJI  THEN 

Il«l  STEP  1  UNTIL  N2  00 

I.I]t>MCNl)^PERTURB*JJ| 

RATION(N2.W.VECTOR.ITCOUNT»INTCH)I 
■1    STEP    i    UNTIL    N2    DO 

■1.2    DO 
tl»Jj3l«VEcT0RCl»l), 

Search  n  columns  of  smtx  corresponoiNg  to  n  ligenvalues 

WITH  POSITIVE  REAL  PARTS.  SMTxl  AND  SMTx2  ARt  ITS  UPPER 

and  lower  halves*  respectivelyi 
I 

STEP  1  UNTIL  N2  DO 
RtU  GTR  0  THEN 

X,«INDEX+1, 
ChCINOEXJIbII 


00028600? 
00028700? 
00028800? 
00028900? 
00029000? 
00029100? 
00029200? 
00029300? 
00029400? 
00029500? 
000296O0? 
00029700? 
00029ft00? 

00029900? 
00030000? 
00030100? 
00030200? 
OOO3O3OO? 
00030400? 
00030500? 
00030600? 
00030700? 
00030800? 
00030900? 
00031000? 
0003UOO? 
00031200? 
00031300? 
00031400? 
00031500? 
00031600? 
00031700? 
00031000? 
00031900? 
00032000? 
00032100? 
00032200? 
00032300? 
00032400? 
00032500? 
00032600? 
00032700? 
00032800? 
00032900? 
00033000? 
00033100? 
00033200? 
00033300? 
00033400? 
00033500? 
00033600? 
00033700? 
00033*00? 
00033900? 
00034000? 
00034100? 
00034200? 
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FOR  J I «i 

STEP  1 

UNTIL 

N  00 

FOR  Ii«l 

STEP  1 

UNTIL 

N  DO 

BEGIN 

SMTXltl 

, J) l-SMTXf I.SE, 

SliTX2(I 

,JJ|»SMTX(N*I 

•SEARCHIJIJJ 

endi 

INVtR$E(N 

.SMTX1, 

>TEMP)I 

MULT(N*N» 

N»SMTX2#TEMP» 

K), 

END  RlCCATl2» 

00034300? 
00034400? 

00034500? 
00034600? 
00034700? 
00034800? 
00034900? 
00035000? 
00035100? 
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APPENDIX  II 


ILLIAC  IV  PROGRAMS 
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*  two  iumc  iv  glypnir  programs.  »clinsys"  ano  "invitation"  are 

I  LISTED  BELOW,  "CLlNSYS"  IS  USEO  TO  OBTAIN  THE  INVERS*  OF  A  MATRIX 

x  ano  the  solution  or  *  complex  system  of  unsymmetric  linear  eouations. 
%   "INVITERATION"  is  useo  to  calculate  the  principal  vectors  of  degree 

%    K  FKOM  (<WLAMBDA*I)**K)*X«0. 

SUBROUTINE  cLINSyS(CINT  N.CINT  R#PcPOINT  A'PCPOINT  x»PCPOINt  B)| 

begin  |  the  subroutine  solves  a  complex  system  of  unsymmethic  linear 
%   Equations  ax«r.  a  is  n*2n  and  b  is  n*?r  matrices*  respectively, 

*  THElR  gEnErAi  EiEMENTS  ARE  All»2J-U*I*ACl»2j]  and  BCI«2J"11* 
X  I*B[I,ZJ],  RESPECTIVELY,  R  IS  THE  NUMBER  OF  RI<*HT-MAND  SIDES. 

SUBROUTINE  INNPROOf CREAL  CNCREAl  C2»PREAL  AK.PREAL  «K#CREAL  01. 
..  CREAL  02)» 

beg1n  i   acCUmulatEs  the  sum  of  products  and  adds  it  to  the  initial 
x  value  (c1.c2), 

CHEAL  Sl2| 

S12!»cUC2| 

S12I»S12*R0*'SUM<AK*B»<)I 

DH.S12, 

D2I.S12-01 I 
END* 

SUBKqUTINE  CDEC0MP0SE(ClNT  N*PCPnINT  A'CREAL  DETR'CRtAL  UEtI* 

CINT  DETE#CNP0INT  INT) I 
BEGiN  %    THE  CUMPLEy  UNSyMMETRIC  mATRIx  A  IS  DECOMPOSED  INTO  LU. 

%   where  l  is  lower  triangular  and  u  is  unit  upper  triangular 
%   matrices,  and  overwritten  on  a.  a  is  stored  in  h*2u   array 

X  And  ITS  GENERAL  ELEMENTS  ARE  At  I»2j*l)*I*A(  I»2j3.  INTIH 
*  K^EpS  ThE  RECORD  Or  ANY  INTERCHANGES  MApE  To  THE  ROWS  Of 
A.  THE  DETERMINANT  ( DETR* I #0ET I  )*2**0ETE  IS  ALSO  COMPUTED. 

K.L.P.PPj 

TOR  ATt64ll 


C1N 
PHE 
PKE 
CHE 
LAB 
MUD 
MUD 

Luu 

I 
DtT 
LUO 

RLG 
L 

L 
B 


T  I»J» 
AL  V£C 
AL  Z7I 
AL  X.Y 

EL  FAI 
EI-TRU 
E«-TRU 
P  T  lal 

NNpRoD 
R  I  •  1  | 
P  Kl.l 
IN 

l»Kl  P 

OOP  1 1 

EGIN 

MpDEI 

MoDE« 

iNNpR 

INNPR 
INNpR 
INNpR 
Y|«-H 

MODE! 
MOOEI 
MDDE« 


.Z. V.W.H.HHJ 

Ll 

El 

E  ANO  PEN  LSS  N*NI 

•1»N  DO 

(0»0»  At  I 

DETII-Oj 

»1»N  00 


J»ACI  J.  INTII  ].W)  I 
DETEI»0| 


I»K*KJ  p 
■K»1.N  D 


PI«P-1I  ZI'OI 

0 


■TRUEl 

■TRUE  AN 

0DCGRA8 

00(-H,-h 

ODCGRAp 

ODCH.HH, 

J 

•TRUEl 

■TRUE  ANO  PEN  EQL  PP*ll   AtI]l«X| 

■TRUEl 


D  PEN  LSS  N*N| 

ONE(AtI].PP-l).O.Ain.AT(PPJ.H»HH^ 
H»RTL(W.Atl])»ATtP].X,HH)l 

ONE(ACn.P-l)»0»RTL(l»»AtI]).ATtPP)»H»HH)l 
All J.ATtp]»H»HH)l 


oooooioo? 

00000200? 
00000300? 
OOOOOaOO? 
00000500? 
00000600? 
00000700? 
00000800? 
00000900? 
00001000? 

ooooiioo? 

00001200? 
00001300? 

oooouoo? 

00001500? 

oooouoo? 
00001700? 

00001800? 
00001900? 
00002000? 
00002100? 
00002200? 
00002300? 
00002400? 
00002500? 
00002600? 
00002700? 
00002800? 
00002900? 
0000300-''? 
00003100? 
000Q32C0? 
C000330A? 
00003^00? 
000035C0? 
00003600? 
00003700? 
00003800? 
00003900? 
00004000? 
00004100? 
00004200? 
00004300? 
00004400? 
00004500? 
00004600? 
00004700? 
00004800? 
00004900? 
00005000? 

00005100? 
00005200? 
00005300? 
00005400? 
00005500? 
00005600? 
00005700? 
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Lll 


M00EI»TRUE  AND  PEN  EOL  P-l J   AC  1 3  a*Yl 

MODEl«TRUE| 

modei«true  AND  PEN»II 
X|«(X*X*y*V)/lNTlIJ| 

If   y  gtR  z  then 

BpGjN  ZI«X|  li*\>    ENDl 
END) 

M00EI,TRUE|  MOOEI.TRUE  AND  PEN  LSS  N+Nl 
IF  L  NEO  K  THEN 
BEGIN  OETRI--OFTRJ 

OETII.-DFTU 

22  a  »At  K  5 1 

AfK]  I-aC.  H 

AtLJI»ZZ| 

INT[L] l«!NTcKj| 
END| 

INT[K]|«LI 

Xl»GRABONE(AtK],pp-l)| 

Yl«GRABONE(A[K],P-i)| 

Z|»X*X*Y*Y| 

Ht«X*DETR-Y*DETI» 

DETf t»X*DETI*Y*DETRl 

DETrIbN| 

IT    ABS(DETR)    GTR    ABS(DETI) 

TmEn   wi»oetr  else   WI«dETII 
IF    w    .0    THEN    BEGIN    DETE|«0| 
BEGIN    LAREL    Lll 
IF    ABS(W)    GEQ    1    THEN 
BEGIN    W|.W*0. 06251 

DETll"DETI*0t0625l 
DETEI»DETE*«» 
GO   TO    Lll 
END| 
ENDl 

BEGIN    LABEL    L2l 
IF    ABs(W)    LSS    0,0625 
BEGIN    WI.WM6I 

DETRl.DETR*16l 
DETII»DETI*U> 
DETEl-nETE**l 
GO    TO    L?» 
ENOi 
ENDl 
LOOP 
BEGIN 

P|»J*J|    PP|«p. 
MnDEUTRUEl 

MflOEl.TRUE    AND    PEN 
lNNpROD( 


GO    TO    FAILI    ENOl 


THEN 


JI«K*l.l»N    DO 
II 


LSS    K*K-2| 


iNNpR0D(-GRA80NE(AtKJ.PP-l).0.AtK],ATtPPl.H.HHJl 
INNPRnD(.H.-HH.RTL(i..ACK]),AT[p;,v,Hi), 

NlJpRbD(-6RABONE(ACKl.P-n;oJlTL    1.    AtPPn.ATrPPi    M    uun 
iNNpROOCH.HH.AtKJ.ATCPJ.H.HH)^  C        ]       'HH)' 


W|«-Hl 

MoDEI«TRUEl 

MnOEl"TRUE  AND  PEN«PP 

A[K],»(V*X*W*YJ/Z| 
MnOEl-TRUEl 


II 


00005800? 

00005900? 

0000^000? 

00006100? 

0000620G? 

00006300? 

00006400? 

000065C0? 

0000660C? 

0000*700? 

00006800? 

00006900? 

0000^000? 

00007100? 

00007200? 

00007300? 

00007400? 

00007500? 

00007600? 

00007700? 

00007800? 

00007900? 

00008000? 

00008100? 

00008200? 

OOOO83OO? 

OOOO84OO? 

00008500? 

00008600? 

00008700? 

00008800? 

00008900? 

00009000? 

00009100? 

00009200? 

00009300? 

00009400? 

00009500? 

00009600? 

00009700? 

0000980O? 

00009900? 

00010000? 

00010100? 
00010200? 
00010300? 
00010400? 
00010500? 
00010600? 

00010700? 
00010800? 
00010900? 

00011000? 

O0O1H00? 
00011200? 
00011300? 
00011400? 
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MoOtl«TRUE  AND  PEN»P-1» 
AfK)l»(K*X-V*Y)/ZI 

END| 
ENO| 
FAILI 
END)  «Or  CDECQMPOSE. 

SUBKOUTINE  CACCSOLVECCINT  N.CINT  R»PCPOINT  A.PCPOINT  AA.tsiPOlNT 

PCPOINT  B»PCP0INT  X#PCPOINT  B8»ClNT  L) % 
BEfiiN  *  SOLVES  AX«B»  *hErE  A  IS  AN  N*2N  COMPLEX  UNSY"METHIc  AN[) 

*  B  IS  AN  N*2R  COMPLEX  MATRICES*  RESPECTIVELY,  AA  J$  LU 

*  DECOMPOSITION  PROOUCEO  BY  "CDECOMPOSE".  THE  KESIUUALS 
%    B8"B-Ax  ARE  ALSO  CALCULATED  AND  AD«B8  IS  SOLVED.  OVER- 

*  WrITINg  0  ON  BB.  i    IS  THE  Num»Er  Qf    ITERATIONS. 


SUBROUTINE  CS0LVE(CINT  N.ClNT  R»PCPOlNT  A.CNPOlNT 

PCPOINT  B)| 
BLGIN 

CINT 
PREA 
CREA 
PREA 

bool 

Loup 

If 


INT , 


I. J.K.KK.P.PPl 
L  VECTOR  BTt6«J» 
L  X.Y#Z»ZUZ2»H.HH| 
L  VECTOR  Ct6«l| 
EAN  AMODEI 

I i«l »i*n  Dn 

iNTtll  NEQ  I  THEN 
LOOP  Jt»R*R.*l»l  DO 

begin 

Xl»GRABONE(BCI],J-l)J 

MODEI-TRUE  ANO  PEN"J-1| 

Bt  13  i«GRABONE(B( INTC11]»J-1)I 

BClNTtUll'XI 

endi 

LOOP    K««R*R»-2.2    DO 
BEGIN 

KkI.K-1  I 

LOOP    I  IM  •  1  »N    DO 

BEGIN 

MODEl-TRUE    AND    PEN    LSS    I* 1-2 1 

lNNPROD(-GRABONE(B(I)*KK-I)»0*AtI)»BTCKK]»H»HH)| 

lNNPR0DCH,-HH»RTL(l#tAtl]).BTCK]»X.HH)l 

INNPROD(-GRA80NE(BCI]»K-1)»0»RTL(1»»A[I1)»BTIKK]'H#HH)I 

lNNPROO(H«HH*A[I)*BT[K]*HtHH)l 

Y|«-HJ 

P|»UII    PP|«P-ll 

Zll«GRABONE( At I)»PP-1)I 

Z2I«GRAB0NE(ACI)»P-1 )l 

Zi«Zl*Zl*Z2*Z2l 

MOOEI«TRUE  AND  PEN.KKI 
BCI  Jl.(X*Zl*Y*Z2)/Zl 

MODEI-TRUE  AND  PEN»KI 
B[I]I«(V#Z1-X*Z2)/ZI 
END; 

LOOP  Il"N.-l,l  DO 
BEGIN 

AMODEI*  TRUE  AND  PEN  GEO  N| 

MnDE«"AMODE» 

lNN»»RGD(-GRABONE(Bt  I  J.KK-p.OtAt  Il»BTtKKJ#H.HH)| 


00011500? 
00011600? 
00011700? 
00011800? 
00011900? 
00012000? 
00012100? 
00012200? 

00012300? 
00012400? 
00012500? 
00012600? 
00012700? 
00012800? 
00012900? 
00013000? 
00013100? 
00013200? 
OOOI33OO? 
00013<»00? 
00013500? 
00013600? 
00013700? 
00013800? 
00013900? 
00014000? 
00014100? 
000U200? 
00014300? 
00014400? 
OO0U50O? 
00014600? 
00014700? 
00014800? 
00014900? 
00015000? 
00015100? 
00015200? 
00015300? 
00015400? 
00015500? 
00015600? 
00015700? 
00015800? 
00015900? 
00016000? 
00016100? 
00016200? 
00016300? 
00016400? 
00016500? 
00016600? 
00016700? 
00016800? 
00016900? 
00017000? 
00017100? 
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lNNPROO(»H,«HH.RTL(l»»A(!l>fBTtK]»H,HH)f 

MoDEl'TRUE  AND  PEN»KK«1I  C  C I  J  I  »H  I 

MQDCl«AHOOei 

lNNpROO(aGRABONE(B[I]*K-l)»0*RTL(l»»ACI])(BTlKK]*H'HH)l 

INNPR00(H.HH»ACI],BTCK](H.HH}I 
MOOEl«TRuE  AND  PEN«K-ll  C[!]I«"HI 
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